library(haven)

###################################################################################################
##### Loading data for study 1 (called "study1_replication_data") and save in object "mydata" #####
##### use the path where you haved saved the data file "study1_replication_data" ##################
##### by uncommenting the line below and provide your path to the file ############################
###################################################################################################

#mydata <- read_dta(".../study1_replication_data.dta")

#################################
##### Recoding and renaming #####
#################################

library(dplyr)

mydata$condition.f <- as.factor(mydata$condition)
mydata$condition.f <- factor(mydata$condition.f,
                    levels = c(1,2,3,4),
                    labels = c("High comp.", "Low comp.", "High warmth", "Low warmth"))

names(mydata)[names(mydata) == "competence"] <- "Competence"

names(mydata)[names(mydata) == "warmth"] <- "Warmth"

#######################################
##### Filtering data: Danish data #####
#######################################

mydata$employment <- as.factor(mydata$employment)
mydata$education <- as.factor(mydata$education)

data_dk <- mydata %>%
dplyr::filter(country==1)

robust_competence_dk <- mydata %>%
dplyr::filter(country==1) %>%
dplyr::filter(competence_true==1)

robust_warmth_dk <- mydata %>%
dplyr::filter(country==1) %>%
dplyr::filter(warmth_true==1)

###################################
##### Filtering data: US data #####
###################################

data_us <- mydata %>%
dplyr::filter(country==0)

robust_competence_us <- mydata %>%
dplyr::filter(country==0) %>%
dplyr::filter(competence_true==1)

robust_warmth_us <- mydata %>%
dplyr::filter(country==0) %>%
dplyr::filter(warmth_true==1)


####################################################################################
##### Figures for main results of study 1 (Figures 2 and 3 in main manuscript) #####
####################################################################################

##### Loading packages needed for visualization #####
library(ggplot2)
library(dotwhisker)
library(cowplot)

##### Panel A of Figure 2: study 1 (DK results for competence) #####
  
##### Calculating point estimates and standard errors via regression #####
dk_admtrust_comp <- lm(admtrust~Competence, data=data_dk)
r_comp_dk_admtrust <- (lm(admtrust~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_dk))

dk_poltrust_comp <- lm(poltrust~Competence, data=data_dk)
r_comp_dk_poltrust <- (lm(poltrust~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_dk))

dk_trustcs_comp <- lm(trustcs~Competence, data=data_dk)
r_comp_dk_trustcs <- (lm(trustcs~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_dk))

dk_trustwas_comp <- lm(trustwas~Competence, data=data_dk)
r_comp_dk_trustwas <- (lm(trustwas~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_dk))


  ##point estimates
  coef.matrix <- matrix(c(0.003, 0.06,
				0.0024, 0.045,
				0.011, 0.11,
				0.02, 0.071), nr=1)
                    
                          
  ##standard errors 
  se.matrix <- matrix(c(0.016, 0.023,
				0.018, 0.024,
				0.02, 0.028,
				0.02, 0.025), nr=1)					
                        
  
  ##variable names
  varnames<- c("Competence (DK respondents)")
  
  model_names <- c("Adm. inst.",
                   "Pol. inst.", "Civil servants", "Gov. in Cop.")

  submodel_names <- c("Exp. manipulations","Attentive subjects")
  
  
  results_df <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix),
                           std.error=as.vector(se.matrix),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))

  
 
  
mainplot1_DK_competence <- small_multiple(results_df) +
    scale_x_discrete(limits = model_names) +  
    theme_bw() + ylab("Marginal effects on citizen trust") +
    geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
    theme(axis.text.x  = element_text(size=12),
	    axis.text.y = element_text(size=10),
	    axis.title.y = element_text(size=12),
          legend.position=c(.99, .99), legend.justification=c(1, 1), 
          legend.title = element_text(size=12),
          legend.background = element_rect(color="black"),
          legend.key.size = unit(12, "pt")) +
	    scale_y_continuous(breaks=c(-0.10, -0.05, 0, 0.05, 0.10, 0.15, 0.20),limits=c(-0.10,0.22)) +
          scale_colour_grey(start=0.2, end=0.8, name="Model") +     
          ggtitle("") 


mainplot1_DK_competence


##### Panel B of Figure 2: study 1 (DK results for warmth) #####

dk_admtrust_warm <- lm(admtrust~Warmth, data=data_dk)
r_warmth_dk_admtrust <- (lm(admtrust~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_dk))

dk_poltrust_warm <- lm(poltrust~Warmth, data=data_dk)
r_warmth_dk_poltrust <-(lm(poltrust~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_dk))

dk_trustcs_warm <- lm(trustcs~Warmth, data=data_dk)
r_warmth_dk_trustcs <- (lm(trustcs~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_dk))

dk_trustwas_warm <- lm(trustwas~Warmth, data=data_dk)
r_warmth_dk_trustwas <- (lm(trustwas~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_dk))

  ##point estimates
  coef.matrix <- matrix(c(-0.014, -0.026,
				-0.009, 0.016,
				-0.007, -0.021,
				0.008, 0.03), nr=1)
                    
                          
  ##standard errors 
  se.matrix <- matrix(c(0.015, 0.027,
				0.018, 0.03,
				0.019, 0.034,
				0.018, 0.03), nr=1)					
                        
  
  ##variable names
  varnames<- c("Warmth (DK respondents)")
  
  model_names <- c("Adm. inst.",
                   "Pol. inst.", "Civil servants", "Gov. in Cop.")

  submodel_names <- c("Exp. manipulations","Attentive subjects")
  
  results_df <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix),
                           std.error=as.vector(se.matrix),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))

  
 
  
mainplot1_DK_warmth <- small_multiple(results_df) +
    scale_x_discrete(limits = model_names) +  
    theme_bw() + ylab("Marginal effects on citizen trust") +
    geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
    theme(axis.text.x  = element_text(size=12),
	    axis.text.y = element_text(size=10),
	    axis.title.y = element_text(size=12),
          legend.position=c(.99, .99), legend.justification=c(1, 1), 
          legend.title = element_text(size=12),
          legend.background = element_rect(color="black"),
          legend.key.size = unit(12, "pt")) +
	    scale_y_continuous(breaks=c(-0.10, -0.05, 0, 0.05, 0.10, 0.15, 0.20),limits=c(-0.10,0.22)) +
          scale_colour_grey(start=0.2, end=0.8, name="Model") +     
          ggtitle("") 


mainplot1_DK_warmth



##### Panel A of Figure 3: study 1 (US results for competence) #####

## Calculating point estimates and standard errors via regression ##
us_admtrust_comp <- lm(admtrust~Competence, data=data_us)
r_comp_us_admtrust <-(lm(admtrust~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_us))

us_poltrust_comp <- lm(poltrust~Competence, data=data_us)
r_comp_us_poltrust <-(lm(poltrust~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_us))

us_trustcs_comp <- lm(trustcs~Competence, data=data_us)
r_comp_us_trustcs <-(lm(trustcs~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_us))

us_trustwas_comp <- lm(trustwas~Competence, data=data_us)
r_comp_us_trustwas <- (lm(trustwas~Competence + age + leftright + polint + soctrust + gender + employment + education, data=robust_competence_us))


##point estimates
 coef.matrix1 <- matrix(c(0.046, 0.099,
				0.039, 0.10,
				0.062, 0.134,
				0.03, 0.06), nr=1)	

##standard error of matrix
  se.matrix1 <- matrix(c(0.023, 0.026,
				0.025, 0.027,
				0.03, 0.035,
				0.02, 0.026), nr=1)


##variable names
  varnames<- c("Competence (US respondents)")
  
  model_names <- c("Adm inst.",
                   "Pol inst.", "Civil servants", "Gov. in Was")

  submodel_names <- c("Exp. manipulations","Attentive subjects")

robust <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix1),
                           std.error=as.vector(se.matrix1),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))


mainplot1_US_competence <- small_multiple(robust) +
    scale_x_discrete(limits = model_names) +  
    theme_bw() + ylab("Marginal effects on citizen trust") +
    geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
    theme(axis.text.x  = element_text(size=12),
 	    axis.text.y = element_text(size=10),
	    axis.title.y = element_text(size=12),
          legend.position=c(.99, .99), legend.justification=c(1, 1), 
          legend.title = element_text(size=12),
          legend.background = element_rect(color="black"),
          legend.key.size = unit(12, "pt")) +
          scale_y_continuous(breaks=c(-0.10, -0.05, 0, 0.05, 0.10, 0.15, 0.20),limits=c(-0.10,0.22)) +
          scale_colour_grey(start=0.2, end=0.8, name="Model") +  theme(plot.title = element_text(hjust = 0.5, face="bold")) +     
          ggtitle("") 

mainplot1_US_competence


##### Panel B of Figure 3: study 1 (US results for warmth) #####

us_admtrust_warm <- lm(admtrust~Warmth, data=data_us)
r_warmth_us_admtrust <- (lm(admtrust~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_us))

us_poltrust_warm <- lm(poltrust~Warmth, data=data_us)
r_warmth_us_poltrust <- (lm(poltrust~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_us))

us_trustcs_warm <- lm(trustcs~Warmth, data=data_us)
r_warmth_us_trustcs <- (lm(trustcs~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_us))

us_trustwas_warm <- lm(trustwas~Warmth, data=data_us)
r_warmth_us_trustwas <- (lm(trustwas~Warmth + age + leftright + polint + soctrust + gender + employment + education, data=robust_warmth_us))

##point estimates
 coef.matrix1 <- matrix(c(0.03, 0.058,
				0.028, 0.066,
				0.038, 0.071,
				0.035, 0.070), nr=1)	

##standard error of matrix
  se.matrix1 <- matrix(c(0.023, 0.027,
				0.024, 0.027,
				0.028, 0.035,
				0.022, 0.025), nr=1)


##variable names
  varnames<- c("Warmth (US respondents)")
  
  model_names <- c("Adm inst.",
                   "Pol inst.", "Civil servants", "Gov. in Was")

  submodel_names <- c("Exp. manipulations","Attentive subjects")

robust <- data.frame(term=rep(varnames[1], times=2),
                           estimate=as.vector(coef.matrix1),
                           std.error=as.vector(se.matrix1),
                           model=as.factor(rep(model_names, each=2)),
                           submodel=rep(rep(submodel_names, each = 1), times = 4))

mainplot1_US_warmth <- small_multiple(robust) +
    scale_x_discrete(limits = model_names) +  
    theme_bw() + ylab("Marginal effects on citizen trust") +
    geom_hline(yintercept = 0, colour = "black", linetype = "dashed") +
    theme(axis.text.x  = element_text(size=12),
 	    axis.text.y = element_text(size=10),
	    axis.title.y = element_text(size=12),
          legend.position=c(.99, .99), legend.justification=c(1, 1), 
          legend.title = element_text(size=12),
          legend.background = element_rect(color="black"),
          legend.key.size = unit(12, "pt")) +
          scale_y_continuous(breaks=c(-0.10, -0.05, 0, 0.05, 0.10, 0.15, 0.20),limits=c(-0.10,0.22)) +
          scale_colour_grey(start=0.2, end=0.8, name="Model") +  theme(plot.title = element_text(hjust = 0.5, face="bold")) +     
          ggtitle("") 

mainplot1_US_warmth


################################################################################
##### Combining panels and plotting Figures 2 and 3 in the main manuscript #####
################################################################################

####################
##### FIGURE 2 #####
####################

fig_2 <- plot_grid(mainplot1_DK_competence, mainplot1_DK_warmth, ncol=1, labels="AUTO")
fig_2

#pdf("plot_DK_study1.pdf")
#fig_2
#dev.off()

####################
##### FIGURE 3 #####
####################

fig_3 <- plot_grid(mainplot1_US_competence, mainplot1_US_warmth, ncol=1, labels="AUTO")
fig_3

#pdf("plot_US_study1.pdf")
#fig_3
#dev.off()
